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ABSTRACT 

Four  GLYPNIR  subroutines  which  may  be  used  in  applications  of  time 
series  analysis  are  the  following: 

(i)   a  subroutine  named  FHANN64  which  computes  the  zero-phase  inpulse 
response  of  a  band-pass  filter, 
(ii)   a  subroutine  named  FFT64  which  computes  the  fast  Fourier  transform 
of  a  real  function, 
(iii)   a  subroutine  named  IFt64  which  computes  the  inverse  fast  Fourier 
transform  (assumed  to  be  real)  of  a  complex  function,  and 
(iv)   a  subroutine  named  AUT064  which  computes  the  autocorrelation  of 
a  series  of  real  values. 

This  report  describes  these  four  subroutines.   Subroutine  declarations  and 
sample  outputs  are  included. 
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INTRODUCTION 

Many  calculations  in  the  area  of  time  series  analysis  are  well  suited 
for  parallel  processing  and  hence  can  be  readily  implemented  on  the  ILLIAC 
IV  computer.   In  particular,  I  have  written  subroutines  which  calculate 
impulse  responses  of  filters  and  fast  Fourier  transforms,  inverse  fast 
Fourier  transforms,  and  autocorrelations  of  one-dimensional  arrays. 
These  subroutines  have  been  written  in  the  higher  level  ILLIAC  IV  language 
called  GLYPNIR  [1,2].   They  will  be  described  and  illustrated  in  the  re- 
mainder of  this  report. 

A  Note  About  Storage.   None  of  the  four  subroutines  to  be  described  in  this 
report  require  the  use  of  matrices.   However,  they  do  use  vector  arrays, 
sometimes  as  input  and  sometimes  as  output  parameters.   In  all  cases  the 
formal  parameters  may  be  either  PE  real  vectors  or,  if  dynamic  storage 
allocation  is  desired,  confined  pointers  used  as  PE  real  vectors.   Further- 
more, all  input  vector  arrays  must  be  in  "straight"  storage  by  rows,  as 
illustrated  in  Figure  1;  all  output  vector  arrays  are  returned  in  straight 
storage  by  rows. 
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FHANN64   Subroutine 

The  FHANN64  subroutine  computes  the  zero-phase  impulse  response  of  a 
band-pass  filter  [3,4,11].   This  is  useful  because  any  signal  may  be 
filtered  in  the  time  domain  through  convolution  with  the  zero-phase  impulse 
response  of  a  given  filter. 


The  band-pass  filter  of  FHANN64  passes  frequencies  f  between  FL  and  FH. 
Its  convergence  properties  are  improved  by  using  Hanning  weights  between 
FLO  and  FL,  and  between  FH  and  FHO  in  the  frequency  spectrum,  which  is 
given  by 


r 


\  (1+  cos    [tt    ^L   _   £>         ]},         0    <   FLO    <   abs    (f)    <   FL 


(FL   -    FLO) 


•  (f)    =    I 


1,    FL    <   abs    (f)    <   FH 


l/i,  r        (f    -    FH)tl 

2    {l+COS    ["    (FHO-FH)]}; 


^_0,    elsewhere. 


FH   <  abs    (f)    <   FHO 


See  Figure   2. 


FLO        FL  FH 

FIGURE   2 


FHO 


>     f 


The  coefficients  of  the  impulse  response  are  given  by 
FILT(m)  =  FHO  +  FH  -  FL  -  FLO 


and 


FILT(m+j)  -  FILT(m-j)  =  h  ( j )  -  w(j),   j  =  1,2,3,...  m 
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where  the  response  is  made  symmetric  about  m  =  (N-l)/2  (n  being  the  number 
of  points  in  the  symmetric  impulse  response).   If  frequencies  are  expressed 
in  Hz  and  the  time  increment  DT  in  millisec,  then 

r 


(j)=< 


Lt/.\    1    (sin    [FHO  x  tt(j)]   +  sin    [FH  x  tt(j)]}   /den,    if  den  4   0 


"[2ttTT]  {cos  [FH  x  tt(j)]}»  if  den  =  °; 


where  DT 

."«>  =  27T  (Iooo)  j 

and  ,    r   /.n    (FHO-FH)  ,2  .  _  . 

den  =  1  -  [tt(j)  x  -i — 1-  ]  j  =  1,2,3, ...,m. 

Similar  expressions,  using  FL  and  FLO  instead  of  FH  and  FHO,  apply  to  w(j). 

The  coefficients  of  the  impulse  response  are  calculated  64  at  a  time, 
one  in  each  PE. 

Use.   The  sine  (subroutine  name  SIN)  and  cosine  (subroutine  name  COS)  system 
subroutines  must  be  declared  prior  to  the  FHANN64  declaration.   The  subroutine 
is  called  as  FHANN64  (N,  L,  FLO,  FL,  FH,  FHO,  DT,  FILT),  where  the  parameters 
are  described  below.   Note  that  by  choosing  suitable  measurement  units, 
spatial  filtering  [11]  in  one  dimension  may  be  done  instead  of  temporal 
filtering. 

The  input  parameters  are  as  follows: 

N  -  the  number  of  points  of  time  in  the  symmetric  inpulse  response. 
N  must  be  an  odd  and  positive  integer. 

L  -  the  maximum  index  of  the  FILT  vector  or  pointer.  There  are  L+l 
rows  of  PE  memory  allocated,  for  the  FILT  vector  array.  L  must 
be  a  non-negative  integer. 

FLO  -  the  low  terminal  frequency,  in  Hz . 

FL  -  the  low  cutoff  frequency,  in  Hz. 

FH  -  the  high  cutoff  frequency,  in  Hz. 

FHO  -  the  high  terminal  frequency,  in  Hz. 

DT  -  the  time  increment  of  the  impulse  response,  in  millisec.   DT 
must  be  positive. 
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The  input  frequencies  must  satisfy  0  £  FLO  <  FL  <  FH  <  FHO.   Also,  we 
require  FHO  x  DT  <  500.   This  ensures  that  FHO  is  not  above  the  Nyquist 
frequency  f  =   0/Tvr/ 7  ,  where  fX7  is  in  Hz  and  DT  is  in  millisec. 


"N    2(DT/     )  '  "N 

The  given  sampling  interval  DT  is  too  long  to  sample  frequencies  above 

V 

The  output  parameter  FILT  is  a  vector  array  consisting  of  the  coeffi- 
cients of  the  impulse  response  in  straight  storage  by  rows,  as  previously 
described  (See  Figure  1) .   It  is  important  to  note  that  the  given  filter 
will  be  added  onto  previous  ones  unless  the  parameter  FILT  is  set  equal  to 
zero  prior  to  the  subroutine  call.   One  can  use  this  additive  feature  to 
form  a  notched  filter,  for  example,  one  which  filters  out  60  Hz  line  noise 
(See  Figure  3).   Example  2  of  the  FHANN64  subroutine  (page  9  )  covers  this 
case. 


10   15 


f(Hz) 


55   60   65 

FIGURE  3 


115   120 


The  FHANN64  subroutine  requires  approximately  5500  words  of  storage. 

The  subroutine  declaration  and  sample  outputs  for  FHANN64  appear  on  the 
following  pages. 
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FFT64  and  IFT64  Subroutines 

The  FFT64  subroutine  is  used  to  find  the  fast  Fourier  transform  [5-12] 
of  a  real  function  in  the  one-dimensional  case.   The  IFT64  subroutine  is 
used  to  find  the  inverse  fast  Fourier  transform  (assumed  to  be  real)  of  a 
complex  function  whose  real  and  imaginary  parts  are  one-dimensional. 
Because  of  the  relatively  high  computational  speed  of  the  fast  Fourier 
transform  for  a  sufficiently  large  number  of  points,  it  may  prove  advan- 
tageous to  filter  a  signal  as  follows: 

(i)   Fourier  transform  the  signal  from  the  time  domain  into  the 
frequency  domain, 
(ii)   Multiply  the  result  by  the  frequency  spectrum  of  the  filter 
(frequency  domain) . 
(iii)   Inverse  Fourier  transform  the  filtered  result  back  into  the 

time  domain.   (Note  that  the  inverse  transform  must  be  real  in 
order  to  correspond  to  a  real  signal  in  the  time  domain.) 
This  procedure  often  takes  less  time  than  performing  a  convolution  in  the 
time  domain. 

The  one-dimensional  Fourier  transform  H(j)  of  a  function  F(k) 

(k=0,  1,  2,  ....,  N-l)  is  given  by 

N-l 
H(j)  =  I     F(k)  exp  (-2irikj/N),  j=0,l, 2, . . . ,N-1, 
k=0 

where  i  is  the  square  root  of  -1.   The  inverse  Fourier  transform  F(j)  of  a 

complex  function  H(k)  (k  =  0,  1,  2,  . ..,N-1)  is  given  by 

1   N-l 

F<J)  =n   I       H0Oexp  (2irikj/N),  j-0,1,2, . . .  ,N-1. 
k=0 

In  either  case  the  Cooley-Tukey  method  [7]  may  be  used  to  compute  the  trans- 
form.  This  method  eliminates  redundant  steps  in  the  calculation  so  that 

2 
only  N  log„  N  complex  operations  are  required  rather  than  N  complex  opera- 
tions.  The  actual  algorithm  used  in  the  FFT64  and  IFT64  subroutines  is 
well  described  in  Reference  6.   However,  before  the  transform  is  performed, 
it  is  required  to  scramble  the  input  array  [6];  this  is  done  in  the  GLYPNIR 
subroutines  by  means  of  an  algorithm  developed  by  Duncan  H.  Lawrie  (see 
listings,  pages  16  and  20)  rather  than  as  described  in  Reference  6. 
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Use  (FFT64) .    All  PE's  must  be  enabled  before  the  subroutine  is  called, 
The  sine  (subroutine  name  SIN),  cosine  (subroutine  name  COS),  and  square 
root  (subroutine  name  SQRT)  system  subroutines  must  be  declared  prior  to 
the  FFT64  subroutine  declaration.   The  subroutine  is  called  as  FFT64  (L,  P, 
F,  HO,  HI,  AMP),  where  the  parameters  are  described  below. 

The  input  parameters  are  as  follows: 

L  -  the  maximum  index  of  the  F,  HO,  HI,  and  AMP  vectors  or  pointers. 
There  are  L+l  rows  of  PE  memory  allocated  for  each  of  these  vector 
arrays.   L+l  must  be  a  non-negative  power  of  two;  so  L  must  be 
0,  1,  3,  7,  15,  etc. 

P  -  the  logarithm  of  (L+l)  to  the  base  two.   P  must  be  a  non-negative 
integer. 

F  -  a  vector  array  comprising  the  real  function  whose  fast  Fourier 
transform  is  to  be  computed. 

If  N  is  the  number  of  points  in  the  Fourier  transform,  then  N  must  be 
a  power  of  two  greater  than  or  equal  to  six  (N  £  64)  for  the  FFT64  sub- 
routine to  work.   The  larger  N  is,  the  more  efficiently  the  subroutine 
operates.   Note  that  we  have  the  following  correspondences: 


L+l 


64 

1 

0 

0 

128 

2 

1 

1 

256 

4 

3 

2 

512 

8 

7 

3 

1024 

16 

15 

4 

etc. 


The  output  arrays  are  as  follows: 

HO  -  a  vector  array  comprising  the  real  part  of  the  fast  Fourier 
transform. 
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HI  -  a  vector  array  comprising  the  imaginary  part  of  the  fast  Fourier 
transform. 

AMP  -  a  vector  array  comprising  the  amplitude  of  the  fast  Fourier 
transform. 

Note  that   F,  HO,  HI,  and  AMP  are  all  in  straight  storage  by  rows  (See 
Figure  1). 

The  FFT64  subroutine  requires  approximately  6500  words  of  storage. 

Use  (IFT64).   All  PE's  must  be  enabled  before  the  subroutine  is  called. 
The  sine  (subroutine  name  SIN)  and  cosine  (subroutine  name  COS)  system  sub- 
routines must  be  declared  prior  to  the  IFT64  subroutine  declaration.   The 
subroutine  is  called  as  IFT64  (L,  P,  HO,  HI,  F) .   The  input  parameters  L 
and  P  have  the  same  meaning  and  usage  as  the  corresponding  FFT64  parameters. 
The  input  parameters  HO  and  HI  are  vector  arrays  comprising  the  real  and 
imaginary  parts  respectively  of  the  function  whose  inverse  Fourier  transform 
is  to  be  computed.   The  output  parameter  F  is  the  vector  array  comprising 
the  inverse  Fourier  transform,  which  is  assumed  to  be  real.   HO,  HI,  and  F 
are  in  straight  storage  by  rows  (See  Figure  1) . 

The  IFT64  subroutine  requires  approximately  6500  +  64  (L+l)  words  of 
storage. 

The  subroutine  delcarations  and  sample  outputs  for  FFT64  and  IFT64 
appear  on  the  following  pages. 
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thts  file  printed  oni  11/16/71  at:  t6:o6»57. 


%%*%%ii%%%%%%%  FAST  FOURIER  TRANSFORM  SUBROUTINE    %%%%%%%%%% 

SUBROUTINE  FFT64(clNT  L»  ClNT  P»  PCPOInT  F.  PCPOInT  HO. 
PCPOINT  HI,  PC  POINT  AMP)j 
QFGIN 
PINT 
C  T  NT 
PRLA|_ 
CRLAl 


OrOW.SROW  J 

I, J .TAB. TAP?; 

G0»G1 »  GOT.GI  T. WO» Wl  »  X; 

pt; 


SCNAMqLE  THE  ImPUT  ARRAY  F. 


*  this  routine  takes  a  p  hit  held  right  iustififo  in  i.  an 

*  inverts  the  ornfh  of  the  hits  therfin. 

PE  INTEGER  SUBROUTINE  TORN  AS  RGA(CINT  P»  PINT  I); 
HtGlN 

IF  P>0  THEN 
HEgIN 
Rr,A:  =  T!  CARlSaPj 

CODF  BEGIN 

L  I  T  (  ?  )    =  1  .  0  »  1  : 

S  h  A  a  r     O(i); 

CSHL(l)    2.il 

CLRA; 

CUR(l)    *C^J 
SHA-tL      i.; 
R  T  A  R        ?', 
TXEFAM(t)  ,  -  3  » 

C  S  H  r  c  n    2a; 
R  T  A  i.      1(1); 
tNf)  CoOE; 
LNn 

ELSF  RGA»=0; 
ENQJ  *  oF  TURN -  —  —  —  —  —  —  •  —  —  - 

T  A  8  2  '  s  P  5 

IF  TAr?>6  THEN  TAh?i=6; 
TAc*  jso-6  ; 

if  tap<o  t^en  ta«:so;   %   thIs  will  cause  the  i.uur  un  y  t 

%     S  T  U  P     AT     1=0, 

Loup   t»=o»i.rev/l(Tab.1)~i    oh 

R  t  b  I  N 

L)rUwJ-SHIFTL(TaR»TuRN(TAR2#PEN)); 

S  R  U  W  J  8 0  R  0 m  +  T  U P  m ( T  A  R  ,  I  )  • 
UR()W«=!)Rf)W+I  ; 

Tab?  ispj 

LOOP     Jl=0.1.63     OU 
i  t  fl  I  N 

IF     TURN(6»SHIFT«(TAK2»SHlFTL(6.nR!H)+PFN))a(PEN+J),t5Mt6] 
THEN    M0[0RUw];apTLfJ»lRiE»FI  SRf"M  ]  ); 
URUw  »=RTL(  1  .MUnF ,  ORllm  ) ; 
S  R  f)  m  l  -  R  T  R  (  1  .  M  0  n  F  .  S  >U)  s)  '. 
LNn  \ 
E'ML); 


x«xxs 

FFT 

1 

%%%%% 

FFT 

73  0 

% 

FFT 

5009 

% 

FFT 

51  OU 

% 

FFT 

5  200 

% 

FFT 

5  3  0  0 

% 

FFT 

5  350 

% 

FFT 

5  'I  0  0 

% 

FFT 

5/450 

% 

FFT 

5500 

% 

FFT 

5550 

% 

FFT 

56  00 

%%%%% 

FFT 

5  7  0  0 

D  %% 

FFT 

5730 

%% 

FFT 

5  7  40 

%%%%% 

FFT 

5  75  0 

%% 

FFT 

5754 

%% 

FFT 

575* 

%% 

FFT 

5760 

%% 

FFT 

57  7o 

%% 

FFT 

5  7  80 

%% 

FFT 

579^ 

%% 

FFT 

5«  oo 

%% 

FFT 

5  8  1  0 

%% 

FFT 

5820 

%% 

FFf 

5  H  30 

%% 

FFf 

5  ft  4  J 

%% 

FFT 

5850 

%% 

FFT 

5«  60 

%% 

FFT 

587  a 

%% 

FFT 

5ft  8  0 

%% 

FFI 

5  ft  9  :i 

%% 

FFT 

5900 

%% 

FFT 

5  9  10 

%% 

FFT 
FFT 
FFT 

59  ?0 
5  9  3,) 
61a^ 

% 

% 

FFT 

6150 

% 

FFT 

6160 

0    < 

FFT 

6170 

% 

FFT 

618  ) 

f, 

FFT 

6  1  9  0 

% 

FFI 

6  2  0') 

% 

FFT 

6210 

% 

FFT 

6??') 

* 

FFT 

6  2  3.0 

* 

FFT 

6?  35 

* 

FFT 

6240 

% 

FFT 

6?50 

% 

FFT 

6?60 

% 

FF  r 

6,370 

t, 

FFT 

6?  80 

I 

FFT 

6290 

* 

FFI 

6  300 

* 

FFI 

6  3  05 
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END  Or  SCRAMBLINb  OF  THE  INPUT  ARRAY  F 

LOUP  T  j  =0  » 1 . L  HO 

BFGlN 

UO*sHo[  I  1 > 

FIRST  ITERATION: 

Q  0  T  :  ■  r,  0  J 

IF  pENt [63:l]=o  THEN 

G0»=GoT+RTL( 1 ,MnUE.GOT ) 

ELSE 
(iOJ«RTR(l»N'ODE,r,OT)-GOTs 

SECONo  ITERATION! 

W  o  I  ■  o ,  0  | 

w  i  » =  - 1  .  0  t 

if   pen. u^n  1  =  0    THEN 

q  E  G  I  N 

«o  i  * l .  r>  i 

w  i :  =  o .  o  J 

FNOj 

got  i  =r,o; 

IF    PEN.  [62:11  =  0     I'HEN 

BEGIN 

(iO»"GOT  +  RTU(Z.MnUE,«jOT)*WOj 

Gi i,RtL(2»M0DE.g0T  )*wl j 

E'NiJ 

ELSE 
*  E G  J  iM 

1,0  :  =  KjP(2»  MHMF.  r,OT  >-GOT*  wO  \ 
GT»a"r,OT*Wll 

ThIRL)  ITERATION: 

PT  s  =  3,1 alb926S?S^9793; 
X:=Pl*pEN» [62J?3*0.?b; 
wo JtCoS(X)  ; 

Wl  ts-sTN(X)  J 

GOT  :  =  r,  0  J 

g  i  t  %  •  r,  1  j 

If  PEN.  [61  :  t  |  =  o  THEN 
q  t  G  J  N 

Go:sGoT  +  RTL(^.MnL)E»bOT->*^o-KrL(a,MnoE.un)*wi; 

(j i  :  =  GlT+RTL(y4.  MODE. GOT)*  s'i  +  RTL*  <*»MOOE.r,l  I  )*wOJ 

END 

ELSE 
3 1  G  J  N 

(i6l«KT-RC*»M00F.#r,0T)-G0T  +  W0  +  GlT*Wll 

(ii:  =  RTR(/»»MnUE.r,lT)-GOT*^l-olT*wu: 

ENl>| 

EnURTn  ITERATION: 


X 

FFT 

% 

FFT 

* 

FFT 

* 

FFT 

X 

FFr 

* 

FFT 

* 

FFT 

X 

FFT 

X 

FFT 

% 

FFT 

* 

FFt 

X 

FFT 

I 

FFT 

% 

FFT 

% 

FFT 

% 

FFT 

* 

FFT 

% 

FFT 

% 

FFT 

t, 

FFT 

% 

FFT 

% 

FFT 

% 

FFr 

I 

FFT 

I 

FFT 

% 

FFT 

% 

FFT 

X 

FFT 

% 

FFT 

% 

FFT 

% 

FF  T 

% 

FFT 

% 

FFT 

% 

FFT 

% 

FFT 

% 

FF^ 

% 

FFT 

% 

FFT 

% 

FFT 

* 

FF  i 

% 

FFT 

2 

FFT 

% 

FFT 

t 

FFT 

% 

FFT 

% 

FFT 

I 

FFT 

% 

FFT 

% 

FF  r 

% 

FFT 

% 

FFT 

X 

FFT 

% 

FFT 

% 

FFT 

i 

,  FFT 

t 

,  FFT 

t 

,  FFT 
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Xj 

*0 

Wl 

60 
61 
IF 

60 
61 

END 

^fc-Vi  I 

60 

^1 

E^J 


aPUPENi  t61  t  -n*U.  125; 
isCnS(X)i 

I«-SIN(-X)| 

Ti»goj 

Tt=Gt  I 

PtM,[60:l]so    THE1* 
N 

j  =  GnT  +  RTL(H»MnOE.G0T)*W0-KTLrH»MnnE»r,l7)*wl; 
»aGlT*RTL(8.MnuEt(;0T)*«l  +  rtTL(fl»rtnnE»r,lT)*W0i 

ELSE 
N 

l«HTR(8»MnDE»r,0T)-ti0T*W0+iilT*Wl; 
lsRTR(8»M00E.r,lT)-G0T*Wl-iilT*rt0: 


M  F  T  H     ITERATION 


«0 
wi 

60 
61 
IE 

qtttl 
^0 
61 

EMU 

at  6  I 

60 

61 

E'^lJ: 


=  PI*PEN.  [60|fl]*0.0625  5 

iscns(X) ? 

:«-STN(X  ) ; 
T  J  s  r,  0  J 

tisgi  ; 

PE\j. [59*11=0    THEN 
N 

laGnT  +  R  Ti-(  16,m,|I)|.:,G0t)*wO"RtL(  1  6  »  M  (10  f  ,  i»  1  T  >  *  w  1  » 
!s6lT+RTL(16.^UL)E.G0T)*rt1+KTL(16>MGDF»blT)*wOJ 


ELSE 

N 

|=RTR(16»M0Dr»60r )"G0T*WO+61 T*W1 : 

F«GlT)*(iOT*Wl-iilT*W0l 


I  -  n   |   r<  (   1  '-'  »  TUMI 

;  =  ktr( 1 6 »  m  n  n i 


bTXTH   ITERATION): 


x : 

wo 
*1 
60 
61 
IF 

Rttil 
UO 
61 

F  aU 

qLGl 
GO 
61 

END: 

HO 
HI 

Enl)  ; 


=PI*PEN.[59jS]*0.o31?5j 

:  =  C  n  s  c  x ) ; 
» =  -  $  i  n  ( x ) ; 

T  :  s  r,  o ; 
T  j  =  p,  1  | 

pFm,  [58*1  ]so    THE* 
N 

J=6oT+RTL(32t  MUDE»(iOT"i*WO-RTL(32»M0DE»61T  )*Wl  J 
JsGl  T  +  RTI.C  32,  HUDE»  GOT  >  *  v/ 1  +KTL<  3?»  MODF»  61  T  >  *  *i  0  J 


tl.SF 


:sRrR(32,MOOr»60r)-GOT*^0+b1T*rtl; 
:sHTR(32,M00r,(ilTJ-GOT*M,"blT*W()! 

(  I  1  : sGO 5 
[  T 1 1 »  1  1  > 

*   of   Loop 


ITLwAtTONS    SEVFN    THWUijr.H    P  +  6: 


X 

FFT 

13  300 

* 

FFT 

13  0  0  0 

1 

FFT 

13500 

K 

FFT 

13540 

% 

FFT 

13560 

1 

FFT 

13600 

% 

FFT 

13700 

% 

FFT 

1  3  ft  0  0 

% 

FFT 

13900 

% 

FFT 

14000 

I 

FFT 

1  a  1  o  o 

I 

FFT 

14  2  0  0 

* 

FFT 

1430O 

% 

FFf 

1  a  a  o  o 

t 

FFr 

14500 

% 

FFT 

1  a  6  0  J 

% 

FFf 

1  4  7  0  0 

% 

FFT 

1  4H0O 

% 

FFT 

15000 

* 

FFT 

15100 

% 

FFT 

15200 

% 

FFT 

15  2  0  0 

k 

FFT 

15260 

i 

FFT 

15  300 

I 

FFT 

15  4  0  0 

% 

FFT 

15500 

% 

FFT 

15600 

* 

FFT 

15700 

/ 

FFT 

15800 

* 

FFT 

15900 

* 

FFT 

16  0  0  0 

* 

FFr 

16]00 

* 

FFT 

16200 

* 

FFT 

1630O 

t 

FFT 

1  6  a  0  0 

% 

FFT 

1650-) 

% 

FFT 

16  7  0  0 

% 

FFT 

1  6  8  0  'J 

% 

FFT 

16  900 

% 

FFT 

16940 

% 

FFT 

1  6  9  6  0 

* 

FFT 

17000 

* 

FFT 

1710.' 

% 

FFT 

1720^ 

% 

FFT 

17  3  0  0 

ft 

FFT 

17400 

% 

FFT 

175  0') 

i 

FFT 

17600 

l 

FF! 

17700 

* 

FFr 

1 7  ft  o  ;> 

* 

FFr 

17  9  0  0 

% 

FF* 

18000 

i 

FFf 

18  10  0 

i 

FF* 

1A200 

/. 

FFl 

1H240 

l 

FFT 

1  8250 

* 

FF- 

18^60 
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3  EG 


BE 


EM 
ENO 


iE 

IN 

ClN 
PRE 
PIN 
IKJ 
FAC 
LOU 

GIN 
1K2 
LOU 
LOU 

L(i  I  N 
Jt: 
Jl  I 
K|» 
X.a 

WO; 
w  i  ; 

UOT 
GO  J 
GIT 
*j  1  ! 
HOL 
Hit 
HOL 
Hit 

NU  J 
Ik: 
FAC 

i)l 


P>0  THEN 


IK.IKP.IT.JI.JIIKi 
FACTOR; 


T 

AL 

T 
T0RJspI*0.01S62b; 

p  tt  1=1  *  l »p  nu 

I=SHIETL( t . TK) J 
P  J  J  s  0  »  I  K  2  » |  00 

P  T  I  s  0  •  1  »  I  *  -  1  n  'J 


sj  + 
K*  = 
PEN 

rAc 
=  cn 
a"s 

»SM 
sHo 

:=m 

=  H1 

Jll 
Jll 
Jll 

J  IT 


T  J 
JI* 

+  SH 

TOR 
S(  X 
INK 
Ot  J 
[JI 
It  J 
C  JI 
i=G 
t=G 
K3  5 
Kl  : 


IK 
IF 

*K 
)  J 
X) 

I  ) 
IK 

I  1 
IK 
OT 
IT 
■G 
-G 


ilu.d; 

» 

» 

i ; 
s 
i ; 

+  G  0  *  w  u  - .,  1  *  rt  1  ; 
+  G  o  +  w  1  «■  g  1  *  w  o » 

0  T  -  r,  0  *  n  0  +  G  1 .  *  n  1  » 
lT-Go*W  1  -G  1  +  -ioi 


%    K=PEM+64*I 


slK?J 

TUR:»FAcTf)R*n«5l 


END: 


calculate,  amp. 

LOOP  T  * -0  » 1 . L  00 

A  M  P  [  I  1  ;  s  S  Q  R  T  (  H  0  r  I  ]  *  H  0  [  T  ]  +  H  U  i  ]  *  H  1  [  T  ]  )  ; 

*  Op  p  F  T  6  4 


% 

FET 

18300 

% 

FFT 

1  a  4  0  0 

% 

FFT 

lasoo 

* 

FFT 

18600 

%    FFT 

18650 

* 

FFT 

18  7  0  0 

t 

FFT 

18800 

< 

FFT 

1  8  90  0 

% 

FFT 

19000 

% 

FFT 

19100 

i, 

FFT 

19  2  0^ 

% 

FFT 

1  9300 

% 

FFT 

19400 

% 

FFT 

19  5  0  0 

% 

FFr 

19  6  0  0 

% 

FFT 

19  700 

% 

FFT 

19800 

% 

FFT 

1  V  9  0 1.1 

% 

FFT 

2  0  0  0  0 

* 

FFT 

2  0  100 

Z 

FFT 

2  0  2  0  0 

% 

FFT 

2  0  3  0  0 

i 

FFT 

2  0  4  0  0 

% 

EFT 

2  0 1, 0 :  t 

% 

ff  r 

2  0600 

% 

FFT 

2  07  00 

* 

FFT 

2  08  00 

* 

FFT 

2  0  9  0  0 

* 

FFT 

2  1  0  0  •) 

% 

FFT 

2110  0 

% 

FFT 

212  0  0 

* 

FFT 

21300 

% 

FFT 

21400 

% 

FFT 

21500 

% 

FF  r 

21600 

% 

FFT 

21?00 

% 

FFT 

21800 

% 

FFT 

21900 
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THTS  ULE  PrINTEO  UN!  11/16/71   ATI  1  6  t  07  t  43, 


%%%%%%%%%%  INVERSE 

SUBROUTINE  IFT 


BEGIN 


PINT  nROW»SRO* 
CINT  T • J»TAR.T 
PREAL  G0.G1.G0 
CRtAL  P I »  N I N  V ; 
PCPOlMT  Flj 


FAST  FOURIER  TRANSFORM  SUHRO 
fa^ClNT  L»  CINT  P,  PCPOINT  HO 

pcpoint  f  ) ; 


UTINE 

%%%%%%%% 
*     PCPOlN 


AB*; 

T.  Gl T. WO. wl »  <; 


* 
i 

% 

%%%%%%%%% 
%  T  Hi 

9!  A  N  V 

PE 
R  L  G  I  N 

If 

rtEr,I 
KQA 

cud 


SC«AMqLE     THE     IMPUT     ARRAYS    H(->    AND    Hi 


Emu 

fc.No 

ENU, 
Flj 


%  t  %  %  %  %  %  %  %  %  % 
s  routine  t 
fhts  the  or 
*  x %%%%% % %% % 

INjFgEH  SUr 

p>o  THEN 

N 

J=T  J  CAR1  Is 

E  BEGIN 

LIT(?) 
SHABR 
CSHLC 1 ) 
CLRAj 
C0R(1  ) 
S  H  A  ^  l 
RTAR 
T  X  E  F  A  M  ( 
CSHR(  1  ) 
RTAl 
C  0  0  E  i 


***  Turn  subroutine  %%%%% 
AKES  a  p  hit  field  right  jUST 
nFK  OF  ThE  ^ITS  THERFIN. 

ROUTINE  TURN  AS  PGAfCTNT  P »  P 


p; 


I  FIFO  IN 

%%%%%%%% 

INT  I  ), 


=  1  »()•  1  ; 
0(1)? 

l ; 
?\ 
1 )  »-3; 

2a; 

111)} 


ELSE  RG 
%  OF  TURN  - 
=GFTPEB(L+1 


A  :  =  0  ; 


%    Fl  IS  THE  IMAGINARY  PART 
<  INVERSE  FrjuRTFR  Transform 


TaB2*sP? 

If  ta«?>6  then 

TAH»sp-6; 

IF  TAq<0  THEN 


TAHi?  S=6  ; 


TAHj=o;    *  THIS  WILL  CAUSE  THE  LOUP 

%     <;TUP  AT  [  =  0. 
VL<  T AH , 1  1-1  OU 


LOOP  f :=0» 1 ,RF 
BEGIN 

UROw J-SHIF TL(T 

SrOin:  =  orow  +  tur 
URUW  *ar)ROW+I  ; 

Tab2:sp; 

LOUP    J«s0»t»63 
BEGliM 

IF     TUpM(6»SHlFT^(TAf-j2»SHlFTLlft.0R0W)  +  PFN))=(PEM+J). 
THt  N 
rt  E  fi  I  N 


AB»  TdRN ( TAd2»  PtlO  )» 
M  (  T  A  ri  •  I )  | 


OU 


%%%%%%%%%%% 
%%%%%%%%%%% 

T  HI,       t 

% 
% 
% 
% 
% 
i 
% 
% 
% 
% 

%%%%%%%%%%% 

I»  ANO   %% 

%% 

%%%%%%%%%%% 

%% 

%% 
%% 
%% 

%% 
%% 

%% 
%% 
%% 

%% 
%% 
%% 

%% 
%% 

%% 
%% 
%% 
%% 

%% 
......... %% 

OF  THF      X 
% 

* 

* 
ON  I  TO  * 
?, 
% 
% 
% 
% 
% 
* 
% 
% 
f5Btftl      % 

% 
% 


IFT 

IFT 

IFT    1 

IFT     1 

IFT    1 

IFT     1 

IFT    1 

IF1     1 

IFT    1 

IFT    1 

IF  T    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

TFT    1 

IFT    1 

IFT    1 

TFT    1 

IFT    1 

IFT     1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFF    1 

IFT    1 

IFT     1 

IFT     1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFt    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT    1 

IFT     I 
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F'CDRUw]  l"RTLC  J»TRUEiHOCSROw])| 
Fl'CORnW]  i"RTL(,N  TRUE»H1  [SROi  J )  I 

tNO! 

UR0w:s9TL(l»M0nE»0R0W); 

SRUW»aRTR(l»HOoF»SHOW)J 

ENUj 
ENUJ 

LnO  Of  scrambling  of   the  Inhut  ARRAYS  ho  anu  hi 

LDUP  I  :  B0» 1 »L  on 
BE<jIN 

Go:  =  Fr I  1  J 
bl»=FH  n; 

FIRST  ITERATION! 

Q  0  T  I  ■  fi  0  I 

b 1 1  *  -  r,  U 

IF    PEm,  L  63  5 1  ]=n    THEN 

Ht(il  \ 

U0:  =  G0T  +  RTL(  1  •  HHt)E»  iiOT  )  ; 

Gl  *«G1  T  +  RTL(  1  »Ml)UE»r,l  T  ) ; 
EHU 

ELSE 
HtblN 

GO  :sHtR(1» mode. got  J -got ; 
b  1  »  =  K  t  R  (  1  •  M  0  D  E  ,  r,  1  T  )  -  G 1  T  I 

SECONo  ITERATION 


WO 

«1 

IF 

^tbl 

WO 

«1 

E'ML) : 

bO 
bl 
IF 

^Ebl 
bO 
bl 

E.JO 

Rtbl 
bO 

i,l 
E  N  0  . 


8=0, 0> 
J  =  1 ,  n ; 
PEN«  t63* 1  ]»0  r Ht N 

N 

:=1.0; 

a=o, o; 

r :  s  r,  0  ; 
Ti-nl  ; 
pEm.[62:1Jso  ThE-»J 

N 

s3b0T  +  RTL(2.  inUE»GOT )*WO-RlL(?»MOnE»Gl  1  )+wl  ; 
JaGt  T  +  RTI.C2,  MHUE.  bOT  )  *  w  1  +  K  TL  C  2  .  MOOE  »  r,  1  T)  +  wO» 


ELSE 


N 


;  =  RTR(2»MrjnE»r,0T)-G0T*W0+GlT*wl: 
l«RTR<2»MnDF.rTlT)-.(iOT*Wl-,,lT*rtOI 


TmIrD  TTERATIQmj 

Pj  **3,  141592h535«9?  ?3; 
X»=PI*PEN. [6?l?]*0.25; 

« o  i  ■  c  n  s  ( x )  i 


X  ] 

FT 

15100 

X  ] 

FT 

15200 

%   ] 

FT 

15300 

X  1 

FT 

1  5  /J  0  0 

X  ] 

FT 

15500 

X  1 

FT 

15*00 

*  ] 

FT 

15700 

*  ] 

FT 

15  8  0  0 

X  ] 

FT 

15  9  0  0 

*  ] 

'FT 

160  0  0 

*  1 

:ft 

1610  0 

*  1 

FT 

16  2  0  0 

X  1 

FT 

16300 

*  ] 

[FT 

16  4  0  0 

X  ] 

:ft 

1  6  <5  0  0 

*  ] 

FT 

16600 

%     1 

'FT 

16  7  0  0 

X  ! 

rFT 

16  8  0  0 

%  1 

FT 

16900 

X  1 

FT 

1  70  0  0 

*  ] 

[FT 

1710  0 

X  I 

[FT 

17200 

X  } 

[FT 

17300 

*  ] 

[FT 

1  7/.t00 

* 

[FT 

17  5  0  0 

X 

[FT 

1  7  6  0  0 

%    ] 

[FT 

17  7  0  0 

X  ] 

[FT 

17800 

i    1 

'FT 

17  9  0  0 

*  1 

[FT 

18  000 

X  1 

[FT 

18  10  0 

X  ] 

'FT 

18  200 

X  ] 

[FT 

18  3  0  0 

X  1 

[FT 

18  4  0  0 

X  1 

[FT 

lft^Ou 

X  ] 

[FT 

18600 

X  i 

[FT 

18700 

{  1 

\FT 

18800 

% 

[FT 

18900 

% 

[FT 

19000 

%    ] 

[FT 

1910  0 

X 

[FT 

19  2  0  0 

X 

[FT 

19300 

X  1 

[FT 

19400 

X  ] 

[FT 

1950U 

* 

[Ft 

19*00 

% 

[FT 

19700 

I 

[FT 

19800 

% 

[FT 

19900 

t    ' 

[FT 

20000 

%    " 

[FT 

20)00 

X 

tFT 

20200 

*  ' 

[FT 

2  0  3  0  0 

* 

[FT 

2  0  4  0  0 

X 

IFF 

2  0500 

% 

[Ft 

20600 

% 

[FT 

20700 
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"1 

GO 

Gi 

IF 

aEGl 

GO 

Gl 

END 

J  EG I 

GO 

Gl 

EiMU  j 


JaSlN(X)  I 

T:  =  r,o; 
T : »  g  1 5 
pEn.  t6i 1 1 j«o   Then 

N 

l«G0T*HTL(4»M0UE»G0T)*W0-!HLf4»M00E»r,lT)*rtl! 
JxGlT+RTL(^»MTUE.GOT)*wl  +  HTL(^#MnnE»r,lT)*^0? 


ELSE 


N 


:=KTR(^»MnOE.GOT)-GOT*WO+GlT*^l; 
JsKtR(^»HODE,  r1lT)-GOT*Wl-.-1lT*rtOI 


fourth  iteration* 


X; 

«0 

N] 

GO 
Gl 
IF 

qtGt 
GO 
Gl 

E'mO 

HtUl 

GO 
Gl 

F'mUj 


■PI*PEN.[61 $31*0. 123; 
Jscns<x); 
JsStn(X) ; 
T  :=go; 
Ti»r,i  5 
dEiu.[6ouj«o   then 

:-G0T  +  RTL(B»wnUE» GOT  )*WU-KTL(B» MODE* GlT)*wl; 
JaGlT+RTLC«»Mf.lUE»r,On*wl  +  HTL(fl»MnnEiGtT)*wOl 


ELSE 


r\| 


:  =HtR(8.  mode,  GOT  )-GOT*W0  +  GlT*«il  J 
isKTR(8»MODEtGlT)-GOT*wl-blT**0; 


MFTH     ITERATION!: 


Xt 

WO 

w1 
GO 
Gl 
IF 

ntGl 
^0 
Gl 

FN  I) 

■!LGT 
GO 
G1 

FND: 


sPl*PF_Nt[6(.)|/i]*0.0625; 

JaCos(X); 
« *  s  t  m  (  x ) ; 
Tj=goj 

T|=r,l  J 

PEm. [59:11=0    THEM 
N 


:=G0T  +  RTL(l6,M;|0t:.G0l)*^0"KTL<16#Mri0F.GlT)*Wi; 
laGlT*RTL(l6,MUUE»GOT)*w1+HTl.(l6.MO0E»blT)*wOI 

ELSE 
N 

|sHTR(16»Ma0F#G0T)-G0T*rt0,«,GlT*Wl; 
lsHTR( 16»MOnF«GlT)"GOT*Wl"GlT*W0! 


SIXTH    ITERATION: 

X|«PI*PEN.[59iS]*0.03125j 
wo»=cns(x) ; 
W 1 1 «  S  T  N  ( x ) ; 

G  0  T  J  =  G  0  I 
G  1  T  J  »  r,  1  J 


% 

IFT 

20800 

% 

IFT 

20900 

% 

IFT 

21000 

% 

IF: 

21100 

% 

IFT 

21200 

% 

IFT 

21300 

% 

IFT 

21400 

% 

IFt 

215  0  0 

% 

IFT 

216  0  0 

% 

IFT 

21  700 

% 

IFT 

2  1 «  0  0 

% 

IFT 

21900 

% 

IFT 

22000 

% 

IFT 

22100 

% 

in 

22200 

I 

IFT 

22100 

% 

IFT 

22aoo 

% 

TFT 

22500 

t, 

IFT 

22*00 

% 

IFT 
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AUT064  Subroutine 

The  AUT064  subroutine  is  used  to  find  the  autocorrelation  of  a  one- 
dimensional  array  [4,  13,  14].   The  autocorrelation  is  the  time-domain 
analogue  of  the  power  spectrum  of  a  signal.   For  an  array  X(j)  (j=0,l,2, . . .n-1) 

the  autocorrelation  is  given  by 
n-l-i 
C(i)  =   I       X  (j+i)  x  X(j),  i  =  0,1,2,.  ...,n*  -  1. 

d-o 

It  is  required  that  the  number  of  points  in  the  autocorrelation,  n',  be 
less  than  or  equal  to  the  number  of  points  in  the  array,  n.    Also,  for 
negative  lags  C(-i)  =  C(i). 

The  algorithm  used  in  the  GLYPNIR  subroutine  is  essentially  the  same 
as  that  given  in  Reference  13;  however,  AUT064  does  not  access  the  ILLIAC  IV 
registers  directly. 

Use.   All  PE's  must  be  enabled  before  the  subroutine  is  called.   The  sub- 
routine is  called  as  AUT064  (M,  L,  X,  C) ,  where  the  parameters  are  described 
below. 

The  input  parameters  are  as  follows: 

M  -  the  maximum  index  of  the  X  vector  or  pointer.   There  are  M+l  rows 
of  PE  memory  allocated  for  the  X  vector  array.   M  must  be  a  non- 
negative  integer. 

L  -  the  maximum  index  of  the  C  vector  or  pointer.   There  are  L+l  rows  of 
PE  memory  allocated  for  the  C  vector  array.   L  must  be  a  non-negative 
integer. 

X  -  the  one-dimensional  array  whose  autocorrelation  is  desired. 

L  must  be  less  than  or  equal  to  M. 

The  output  parameter  C  is  a  vector  array  consisting  of  the  autocorrela- 
tion coefficients  for  non-negative  lags.   Both  X  and  C  are  in  straight 
storage  by  rows  (See  Figure  1) . 

The  AUT064  subroutine  requires  approximately  5000  words  of  storage. 

The  subroutine  declaration  and  sample  outputs  for  AUT064  appear  on  the 
following  pages. 
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%%t%%tt1>%%%%%%%%%%  AUTOCORRELATION     SOBROOTINt  %%%%%%%%%% 

%%f %%*%%%%%%%%%%%%%%%%%%%%%%%%%%%%*%%%%%%%%%%%%%%%%%%%%%%%%% 

SlJdKnUTlNE    AUTU64rClNT    M,     ClNT    L,    PCPGINT    X,    PCPUINT    C)l 
BEG  I'M 

PRLAL 

CInT 

BOOLEAN 
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XROW. rRUW. TEMPj 
I »  J »  K  •  M I » I J 1 ; 
SA\/EMnot ; 
LO.L^ j 
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ENDj 
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ENDj 

M0DE:=TRUF j 
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IF  L  =  M  THEN  (iO  TU  LM» 

LnnP  ji«o»i.mi-i  on 

R  E  G  I  N 

I Jl  JsL  +  J+1  ; 

LnOP  K»sO» 1 ,63  00 

b  E  G  I  N 

CpnWI»GRABnNt(X[ J]  ,K)*XHUW  +  CWHWJ 

XRnw;aRTL(  1  »MOl)E»XROW)  ; 

IF    PEN=63    TMEN     XRfiw:  =  RTKC63-K.TPUE»  XC  Ul]  )  , 

EMDJ 
EUD« 
LMi       LnnP  KJsO.l  ,63  OU 

SavEMUOE t=HnUE; 
MnnE!sTRijE: 

TrMPlaftRAB  in|E(  At  «I  1  »K  )  I 
MnnE  jsSAvE'inut ; 

CRnWj  =  TEMH*)(BU'N     +     CROrt; 

Xpnw  j  =  rtl(  1 » mine  »xrl)W  ) : 

HfiDtloSHlFTLt  1  »  HOOF)  J 
END} 

MnnE:  =  TRUE : 
c  r  L  1  » s  c  R  n  W  ; 
END?  *  l)F  AUT064  --.----.---------------------------- 
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